Effects of global change on snakebite envenoming incidence up to 2050: a modelling assessment

Summary Background Human activities are driving climate, land cover, and population change (global change), and shifting the baseline geographical distribution of snakebite. The interacting effects of global change on snakes and communities at risk of snakebite are poorly understood, limiting capacity to anticipate and manage future changes in snakebite risk. Methods In this modelling study, we projected how global change will affect snakebite envenoming incidence in Sri Lanka, as a model system that has a high incidence of snakebite. We used the shared socioeconomic pathway (SSP) scenario analysis framework to integrate forecasts across the domains of: climate change (historical trend from WorldClim plus three underlying regional circulation models [RCMs] in the Coordinated Regional Downscaling Experiment-South Asia repository, with two emissions pathways [representative concentration pathways RCP4.5 and RCP8.5]); land cover change (Dyna-CLUE model); and human population density change (based on Gridded Population of the World data) from Jan 1, 2010 to Dec 31, 2050. Forecasts were integrated under three different development scenarios: a sustainability pathway (SSP1 and no further emissions), a middle-of-the-road pathway (SSP2 and RCP4.5), and a fossil-fuelled pathway (SSP5 and RCP8.5). For SSP2 and SSP5, we nested three different RCMs (CNRM-CM5, GFDL-CCM3, and MPI-ESM-LR; mean averaged to represent consensus) to account for variability in climate predictions. Data were used as inputs to a mechanistic model that predicted snakebite envenoming incidence based on human–snake contact patterns. Findings From 2010 to 2050, at the national level, envenoming incidence in Sri Lanka was projected to decrease by 12·0–23·0%, depending on the scenario. The rate of decrease in envenoming incidence was higher in SSP5-RCP8.5 than in SSP1 and SSP2-RCP4.5. Change in envenoming incidence was heterogenous across the country. In SSP1, incidence decreased in urban areas expected to have population growth, and with land cover changes towards anthropised classes. In SSP2-RCP4.5 and SSP5-RCP8.5, most areas were projected to have decreases in incidence (SSP5-RCP8.5 showing the largest area with incidence reductions), while areas such as the central highlands and the north of the country showed localised increases. In the model, decreases occurred with human population growth, land use change towards anthropised classes (potentially shifting occupational risk factors), and decreasing abundance of some snake species, potentially due to global warming and reduced climatic and habitat suitability, with displacement of some snake species. Interpretation Snakebite envenoming incidence was projected to decrease overall in the coming decades in Sri Lanka, but with an apparent emerging conflict with sustainability objectives. Therefore, efforts to mitigate snakebite envenoming incidence will need to consider the potential impacts of sustainability interventions, particularly related to climate and land use change and in areas where increases in incidence are projected. In view of global change, neglected tropical diseases and public health issues related to biodiversity, such as snakebite, should be managed collaboratively by both environment and health stakeholders. Funding UK Medical Research Council.

Global change effects on snakebite: an assessment of environmental and health sustainability trade-offs in Sri Lanka

Climate change downscale
We corrected bias using WorldClim v2 1 as observed climate and the RCM 970-2006 as predicted climate.We extracted WorldClim values from a 15 km radius at the centre of each RCM grid cell, and then calculated the difference between minimum and maximum temperatures (T bias = T observed -T predicted ) and the relative difference for precipitation (R bias = R observed /R prediced ) 2 to avoid negative values.Then we averaged the bias factors by month so that we obtained 12 bias correction factors for each month's temperature and rainfall per RCM and RCP.
To estimate absolute climate changes per year, first we corrected bias of the regional circulation model projections of 2006-2050 with the reverse arithmetic, substituting T predicted and R predicted with the RCM values of 2006-2050, so that bias-corrected climate is T corrected = T predicted + T bias , and R corrected = R bias × R predicted .Then temperature and rainfall changes per month (ΔT and ΔR) with respect to observed historical climate (monthly averages of WorldClim 1 ) are: ΔT = T corrected -T observed , ΔR = R corrected /R observed 2 .
Resulting monthly changes were interpolated with thin plate splines from their original resolution of 55 km to the target resolution of 1 km (ΔT 55 km → ΔT 1 km , ΔR 55 km → ΔR 1 km ) 3 , to match WorldClim's 30" dataset 1 .To improve the interpolations over Sri Lanka, we extended the geographical extent by 5° north and west from Sri Lanka including some values over southern India.Once changes were interpolated, we performed the reverse arithmetic used to correct biases, transforming interpolated changes to temperature and rainfall estimates: T predicted, 1 km = T observed, 1 km + ΔT 1 km ; R predicted, 1km = R observed, 1 km × ΔR 1 km .Using relative changes for rainfall did not guarantee that all values were positive, as thin plate splines can still generate negative values.Hence we forced all negative values to be equal to the bottom 1 st percentile of positive values, and to limit the prediction of extreme values we forced all values larger than the 99 th percentile to be equal to that same value for each year, to allow annual variability of extreme rainfall 2 .
To validate downscales we calculated the root mean squared error (RMSE), correlation coefficient and linear regression intercept and slope for the predictions between 1970-2006 with the matching climate from the Climate Research Unit v4.03 (CRU ) 4 .With the downscaled 1970-2006 RCM predictions we extracted and averaged monthly values from a 15 km radius at the centre of each CRU grid cell from 1 km downscales.Then, we calculated the mentioned statistics.The aim of the described tests was to obtain different measures of the size and type of errors of the predictions: RMSE should be close to zero (0); correlation coefficient positive, close to 1 and significantly different from zero (0); linear regression intercept not significantly different from zero (0) and slope as close as possible to 1 and significantly different from zero (0).Given the lack of available RCP 2.6 in CORDEXSA we used current climate to represent a no-further climate change sustainability pathway (see Figs. 1 and 2).

Land Cover Base land cover layer
We used Landsat images to produce the base land cover layer of Sri Lanka.We selected the best available surface reflectance level 2 images (Landsat 4-5 TM and Landsat 7ETM+) for the year 2010 with cloud and haze cover <10% (dry season) [5][6][7] for our classification from https://earthexplorer.usgs.gov/.In the case of Landsat 7ETM+ and regions with low quality images we filled the stripping gaps, finding the missing values in the images closest in time.Once these images were cleaned, we extracted the optical and NDVI bands to classify into discrete classes.
The classification routine used was the Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) and visual interpretation of the clusters Finally, we validated our classification using the ground truth points obtained by visual interpretation of high resolution Google Earth and Google Street view images and expert opinion.We chose 600 random points around Sri Lanka and drew polygons in areas where we could identify homogeneous land cover (images available for year 2010).The overall accuracy of the classification was ~ 95%.

Land cover change simulation modelling
First, we set class transitions rules derived from observed changes between 2010 -2017 as: forest ↔ degraded forest ↔ agriculture → urban, and depending on the location forest ↔ degraded forest ↔ tea → urban.These transition rules, population growth and the area demanded from each land cover class define the different SSPs: In SSP1 demand of urban decreased, and favoured forest recovery; SSP2, agricultural lands grew but mainly at the expense of degraded forest; and in SSP5 urban areas and agriculture grew at the expense of forest and degraded forest.To represent the effect of governance on land use change dynamics, we forced natural protected areas to experience zero deforestation under SSP1 and SSP2, and to be open to deforestation under SSP5.
Second, annual rates of class changes were estimated with the observed changes in the 2010-2017 period and assume that these remain constant over time.Third, we selected a set of dynamic (that change with time, Table 2) and static (Table 3) environmental variables that were used to model the probability of presence of each land cover class with a logistic regression (hereafter called location factors).Fourth, we estimated class demands using the total area covered by each land cover class using the Land Use Harmonization (LUH) data 10,11 .For consistence with the base land cover layer, we merged the corresponding LUH classes to represent the above five classes (Table S5).The LUH data represent land cover as proportion of each class in 0.25 × 0.25° pixels in a time series to year 2100.Therefore the demand corresponds to the sum of %cover × pixel area.As a result this data consists of a table with area values arranged in five columns (land cover classes) and 41 rows (time steps).To allocate the demanded areas for each class, Dyna-CLUE uses the location factors (Tables S2-3) and class transition rules to decide where and how much land of each class will be converted to another.
Finally, we specified spatial restrictions to land cover transitions with a raster layer of distance to the edge of natural protected areas (for SSP1 and 2).Together, land cover demands, location factors and spatial restrictions define the SSP that are represented by land cover predictions.For instance, climatic dynamic factors relate to climate change via the RCP downscales used in projecting the location factors, while human population growth correspond directly to SSPs.
To determine land cover transition rules we analysed the land cover changes between 2010 and 2017 across Sri Lanka.These transition rules take place interactively with the location factors estimated with logistic regressions.These regressions used the presence or absence of each land cover class in 2010 as dependent variable and the environmental variables in Tables S2-3 as explanatory.The presence or absence data for each land cover class were sampled with randomly generated, spatially-independent points.The final logistic model for each land cover class was selected from a model incorporating all the dynamic and static environmental variables and progressively simplifying it until we minimised the Akaike information criterion.The coefficients estimated for each of the final regression models, are those used by Dyna-CLUE as location factors (Table S4).An example of the importance of location factors is that the transition degraded forest → tea only occurs in the central highlands (Fig S2 ), represented by elevation and rainfall variables.

Tree Cover
We allowed a number of tree cover transitions: 5 ↔ 4 ↔ 3 ↔ 2 ↔ 1; 5 →1, 2, 3; 4 → 3, 2, 1; 3 → 2, 1; and reversible 2, 3 → 4; 3, 4 →5.Reversible transitions were only allowed to occur after a period of 20 years, roughly representing the recovery time of vegetation.In scenarios SSP1 and SSP2, tree cover changes were restricted with a raster layer of distance to NPAs.The demand and consequent location factors were calculated with logistic regressions for each tree cover interval comparing tree cover of 2000 with 2010 and with the same dynamic and static location factors as covariables shown in Tables S2-3.
To make tree cover predictions compatible with the original data used as covariates in the snake models, we calculated a correction factor after reclassifying tree cover categories back to proportions, 1 → 0.1, 2 → 0.3, 3 → 0.5, 4 → 0.7, 5→ 0.9.The correction factor (T Corr ) was the average difference between the reclassified projected layers and observed tree cover of 2010 in logit scale (logit(tc) = log(tc/1tc), tc = tree cover): T Corr = E[logit(tc proj ) -logit(tc 2010 )], so as to force tree cover predictions to be bounded by 0 and 1.
Then the logit-corrected tree cover for year i was tc' i = logit(tc i ) -T corr , and corrected tree cover is Tc i = exp(tc' i )/1 + exp(tc' i ).

Tables Table S1.
Correspondence between shared socioeconomic and representative concentration pathways.

Fig. S3 .
Fig. S3.Predicted envenoming incidence trends between SSPs and among RCMs used (top panels).y-axis scale is national average predicted envenoming incidence.Variability of annual rates of change for each scenario and RCM (bottom panels) where panels indicate the scenario, and colour indicates the climate data or RCM used for projection.

Fig. S5 .
Fig. S5.Percentage of change of the country-wide median of the log distance to each species' niche centroid as a function of time grouped by RCM and SSP-RCP and smoothed (top).The current climate plot has been omitted because it is a straight line.Positive relationships between median distance to the centroid and time implicate that there is lower climatic suitability which translates in lower potential abundance (bottom).

Fig S6 .
Fig S6.Base tree cover layer used for generating tree cover predictions.

Fig. S9 .
Fig. S9.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for for SSP2 in four time steps (rows) in RCM CNRM-CM5.Grey colour indicates total absence of suitable conditions.

Fig. S10 .
Fig. S10.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for SSP2 in four time steps (rows) in RCM GFDL-CM3.Grey colour indicates total absence of suitable conditions.

Fig. S11 .
Fig. S11.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for SSP2 in four time steps (rows) in RCM MPI-ESM-LR.Grey colour indicates total absence of suitable conditions.

Fig. S12 .
Fig. S12.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for SSP5 in four time steps (rows) in RCM CNRM-CM5.Grey colour indicates total absence of suitable conditions.

Fig. S13 .
Fig. S13.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for SSP5 in four time steps (rows) in RCM GFDL-CM3.Grey colour indicates total absence of suitable conditions.

Fig. S14 .
Fig. S14.Snake suitability predictions (potential abundance in log 10 scale) corrected for relative abundance for SSP5 in four time steps (rows) in RCM MPI-ESM-LR.Grey colour indicates total absence of suitable conditions.

Table S2 .
Dynamic location factors used to influence transition probabilities.

Table S3 .
Static location factors

Table S4 .
Regression models based on dynamic and static factors for Dyna-CLUE models

Table S5 .
Equivalence between LUH and base line land cover classes

Table S6 .
Average performance statistics of the downscaling method for each regional circulation model compared with observed climate from the Climate Research Unit.

Table S7 .
Annual snakebite envenoming changes estimated under each shared socioeconomic pathway and global circulation model.The percentage change shown is accumulated until year 2050.